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Using mesa array of BiaS^CaCi^Os single crystal was demonstrated recently as a promising route to en- 
hance the radiation power generated by Josephson oscillations in mesas. We study the synchronization in such 
an array via the plasma waves in the base crystal. First, we analyze plasma oscillations inside the base crystal 
generated by the synchronized mesa array and the associated dissipation. We then solve the dynamic equation 
for superconducting phase numerically to find conditions for synchronization and to check the stability of syn- 
chronized state. We find that mesas are synchronized when the cavity resonance of mesas matches with that of 
the base crystal. An optimal configuration of mesa arrays is also obtained. 

PACS numbers: 74.50.+r, 74.25.Gz, 85.25.Cp 



I. INTRODUCTION 

Soon after the discovery of Josephson Effects, it was real- 
ized that Josephson junction can be used to generate electro- 
magnetic waves. When the junction is biased in voltage state 
with voltage V, the two superconducting electrodes have en- 
ergy difference 2eV. The system is similar to two-energy level 
system in atomic physics. When Cooper pairs tunnel from 
the electrode with higher energy to that with lower energy, a 
photon with angular frequency co = 2eV/ft is emitted. The 
frequency can be tuned by voltage and 1 mV corresponds to 
0.483 THz. The radiation power from one junction however 
is weak, of the order of 1 pW. 1 3 Arrays of Josephson junc- 
tions are fabricated to enhance the radiation power 4 8 . Once 
these junctions are synchronized, the total radiation power is 
proportional to the number of junctions squared. 

A stack of Josephson junctions is naturally realized in some 
layered cuprate superconductor 9 , such as Bi2Sr 2 CaCu208 + ,5 
(BSCCO). Because of the large superconducting energy gap 
(60 meV), these build-in intrinsic Josephson junctions (IJJs) 
may have Josephson oscillations with frequencies in the ter- 
ahertz (THz) band. IJJs are packed on nanometer scale, 
much smaller than THz electromagnetic (EM) wavelength, 
and are homogeneous. The THz generator based on IJJs thus 
is promising to fill the THz gap! 10 ! 11 !. L ots f e ff or t has been 
made to excite the coherent THz radiation experimentally in 
the last decade™. On the theoretical side, numerical simula- 
tions and analytical calculations are performed to understand 
the mechanism of radiation.GSE! 

Coherent radiations from a mesa structure of BSCCO in the 
absence external magnetic fields were observed experimen- 
tally in 2007 25 , which renewed the interest in this field. It was 
found that the mesa itself forms a cavity to synchronize the ra- 
diation in different layers, as evidenced from the dependence 
of the radiation frequency / on the lateral size L of the mesa, 
/ = cq/{2L) with Co = cj yfe[. the Josephson plasma velocity 
where e c is the dielectric constant of BSCCO. The cavity res- 
onance mechanism has been confirmed by many independent 
experiments^" 3 ^ and the radiation power is enhanced to about 
30 yuW. 31 A dynamic state with n phase kink was proposed 



to account for the experimental observations^! 3 ^). It was sug- 
gested that the strong in-plane dissipation is responsible for 
the excitation of cavity mode uniform along the c-axis. 34 

From application perspective, the radiation power in the 
present experimental design is still too weak to be practically 
useful. A natural way to enhance the radiation power by using 
thicker mesas has several challenges. First, for a thick mesa it 
becomes difficult to cool the system efficiently. The dissipa- 
tion hence self-heating increases with the volume of the mesa, 
while the heat removal rate remains the same because the heat 
is mainly removed through the substrate. It has already shown 
experimentally that even for a mesa with thickness of ~ 2/j.m, 
central part of the meas is driven to the normal state by the 
severe self-heating! 27 l 28 [ Secondly, it was calculated that for 
a tall mesa a long-range instability destroying the in-phase 
plasma oscillations develops 35 , and only parts of the mesa can 
be synchronized.^ 

To enhance the radiation power while minimizing the self- 




FIG. 1. (color online) Schematic view of multiple mesas atop of 
BSCCO crystal. The mesas are biased independently. The BSCCO 
sample (blue) is sandwiched by gold electrodes (orange). In analyt- 
ical treatment, we consider an array of identical mesas with width L 
and period a. The top surface of the base crystal is free. In simula- 
tions, we consider two identical mesas with width L separated by a 
distance L,„. The mesas are located at the position L away from the 
edges of the base crystal. The top surface of the base crystal is also 
covered by gold electrodes through which the current is extracted. 
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heating, one may use multiple thin mesas on top of the 
same BSCCO single crystal. The multiple mesa structure has 
been fabricated recently and the radiation power is enhanced 
under appro priate conditions as demonstrated in the recent 
experiments! 36 ! 37 ! The mechanism of synchronization among 
mesas is not known. There are two sources of interaction. The 
mesas interact through the radiation fields. They also interact 
through the plasma oscillations in the base crystal. The reso- 
nance damping due to the leaking of radiation from the mesa 
into the base crystal has been considered in Ref. [38]and it was 
demonstrated that this channel gives the main contribution to 
the dissipation. Therefore, synchronization mediated by ra- 
diation fields inside crystal is probably a dominating mecha- 
nism. The present work is devoted to understanding the syn- 
chronization of multiple mesas through plasma oscillations in 
the base crystal and to finding an optimal configuration for 
synchronization. 



II. MODEL 

We consider arrays of identical mesas with the period a 
atop of BSCCO single crystal as schematically shown in Fig. 
[T] Every mesa contains N m junctions and has width L while 
the base crystal contains N c junctions. No external magnetic 
field is applied to the BSCCO. Each mesa is biased indepen- 
dent by a dc current injected from top of the mesa and ex- 
tracted from the sides of the mesa. In this case, the junctions 
in the basal crystal remain zero-voltage state, and the mesas 
are driven into resistive state. The system is assumed to be 
uniform along the y direction and the problem becomes two 
dimensional. The dynamics of the gauge invariant phase dif- 
ference 9„ and magnetic field h„ in the n-th junction are de- 
scribed b\' 3 ^3 
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where A (2) // = fi+\ + /;_] - 2// is the finite difference opera- 
tor. Here the time and coordinate are measured in units of the 
inverse Josephson plasma frequency 1 /a> p and the Josephson 
length Aj = ys correspondingly and the unit of magnetic field 
is i>()/(27ry.s 2 ), where y is the anisotropy factor and s is the 
interlayer spacing. Here u> p = c/(A c yfe^) and <l>o = hc/(2e) is 
the flux quantum. These reduced equations depend on three 
parameters, v c = 47rcr c /(e c a> p ), v a b - Ancr a bl{s c (i) p y 2 ), and 
( = A a bl s, where cr c and o~ a j, are the quasiparticle conductiv- 
ity, and A c and A a b are the London penetration depth along the 
c-axis and ab-plane respectively. The dimensionless electric 
field is given by E Zi „ = d,9„, with E in unit of &uu) p l{2ncs). 

For the mesa with thickness of 1 fim which is much smaller 
than the wave length of THz EM wave in vacuum, there 
is a significant impedance mismatch between the mesa and 
vacuum 44 . Most part of energy is reflected at the edges of 



mesa and cavity resonance is achieved. We can use the bound- 
ary condition that the oscillating magnetic field vanishes at the 
edges. The boundary condition at the edges of the mesas is 
d x 9„ = ±LI M /2, and the boundary condition at the edges of 
the base crystal is d x 9 n = 0, where I p is the bias current in 
the //-th mesa. We assume that the IJJs stack is sandwiched 
by good conductors, such that the tangential current inside the 
conductor is zero, which corresponds to the boundary condi- 
tion d z h(z) — in the continuum limit. 



III. PLASMA OSCILLATIONS AND ASSOCIATED 
DISSIPATION IN THE SYNCHRONIZED STATE 

In this section, we calculate the plasma oscillations and its 
dissipation in the synchronized state assuming that the array 
contains large number of mesas so that it can be treated as 
an infinite system. The time dependence of the phases in 
mesas in resistive state has the form 9„{x, t) = cot + <p„(x) + 
Re[0„(jc)exp(-/wT)] + ifr^ and in the crystal the phases have 
only have small oscillations. Here if/ p accounts for the phase 
shift between synchronized mesas. We consider the case with 
iff M — and leave the more general case for numerical sim- 
ulation in the next Section. We consider voltage range cor- 
responding to the Josephson frequency close to fundamental 
cavity resonance co as to x = (n/L. For d efinit eness, we assume 
that in mesas the kink state is formed 3233 providing strong 
coupling to the cavity resonance meaning that we can use ap- 
proximations <f n {x) as 7rsgn(x) and sin6>„ = Re[/ exp(-i'o»r - 
iip„(x))] as g(;t;)Re[/exp(-/«T)] with g(x) as sgn(x). However, 
a particular shape of the modulation function g(x) has no im- 
portance in further derivations. For isolated mesa on the top 
of bulk crystal this problem was considered in Ref. l38l where 
it was concluded that leaking radiation into crystal provides 
dominating mechanism of resonance damping. 

The amplitudes of phases and magnetic fields obey the fol- 
lowing equations: in mesas for \x - ma\ < L/2, < n < N m , 

(co 2 + iv c cj)9 n + { 2 -^ =ig(x), (3) 
89 

€ 2 V 2 n h„-{l- iVabO)) h„ + ( 1 - iVabCO) -t£=0, (4) 

and in the crystal, for -7V C < n < 0, the first equation has to 
be modified as 



U-l+iv c oj)9„ + £ 2 ^ =0 



dx 



(5) 



Using presentation g(x) = 8m sin (p m x) with p m = 

(2m - l)n/L, we can find solution of these equations as mode 
expansions. In particular, the oscillating magnetic field in 
mesas can be written as 



ipm(gm+a m cos [q,„ (n-N m - 1/2)]) 

U> 2 + iv c (x) - ftp 1 ^ 



cos (p m x) 
(6) 

for n > 0, where q m = q(p m ,ai) with Im(g,„) > is the wave 
vector describing propagation of the plasma wave along c-axis 
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FIG. 2. (color online) The frequency dependence of the decay length 
along odirection, A^, for the fixed in-plane wave vector and repre- 
sentative parameters specified in the plot. 
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FIG. 3. (color online) Comparison of the frequency dependence of 
the oscillating magnetic field at the mesa top for isolated mesa and 
mesa array. The frequency dependence of the oscillating magnetic 
field at the bottom of the base crystal for mesa array is also shown. 
Parameters used in calculations are specified in the plot. 



for the fixed in-plane wave vector and frequency, 



1 - iv ab u 
cos q m = 1 + — ^ — 1 1 ~ 
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In crystal, n < 0, the oscillating magnetic field can be pre- 
sented as Fourier series, 



hf\x) = - ^ H k cos [q k (n+N c + l /2)] exp(rfcc) (8) 



k=2irl/a 

with cos q k = 1 + — —2 — 1 1 
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It is crucial that the synchronized mesa array excites discrete 
set of modes inside the crystal. For the fixed k the frequency 
range to 2 < 1 + i 2 k 2 corresponds to propagating waves along 
the c-axis while the range co 2 > 1 +( 2 k 2 corresponds to evanes- 
cent waves. At the frequency w 2 = 1 + l 2 k 2 the uniform 
plasma mode is excited. The decay length of the plasma mode 
in terms of the number of junctions, Nd(co) = 1 /\m(qk(pS)), 
has a sharp maximum at this frequency, see Fig. [2] For 
a ^ 2L the mode with the wave vector k = 2n/a plays 
the most important role, because frequency of the uniform 
mode a> = { 2 {2nla) 2 is close to the cavity-resonance 

frequency inside the mesa a>] - InjL. 

The unknown coefficients a m and Hk have to be found from 
matching at the interface, h^\x) - h^(x) for n — 0, 1. Tak- 
ing the projection of the equation h ( ^\x) = h^ix) to mode 
m, using 



rL/2 

1 

J-L/2 



dx cos (p,„x) cos kx 



2(-l)> m cos [kL/2] 



il-k 2 



(9) 



we obtain equation expressing a m via Hk 

(gm + a m cos [q m (N m + 1 12)]) 

CD 2 + iv c u) - fip^ 
1 V w r , 1 Pm cos [kL/2] 

= ~ 2, Hk C0S W* ^ + 1/2) J } 2 T2 ■ 

a k=2nl/a L Pm k 

Note that S m (k) satisfy orthogonality conditions 
2 



aL 



2^ SJk)S m >(k) = S„ 
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On the other hand, the inverse Fourier transform of the equa- 

us to express Hk via a m 

iPm (gm + a m cos [q m . (N m - 1/2)]) 



tion h. (x) = h ( "\x) allows us to express Hk via a 



H, 



cos [q k (N c + 3/2) = ) — : 

to 1 + I 



v c co - e 2 p 2 , 



2(-Y) m p m cos [kL/2] 



(10) 



Pi'k 2 

Eliminating Hk, we obtain the linear equations for a m 
(cos [q m (N m - 1/2)] - cos [q m (N m + 1/2)]) a m 

CO OO 

- Jmrn'Om' COS [q m , (N m - 1/2)] = ^ Jmm'gm' (1 1) 



m'=l 

with the matrix 



m' = l 



p,„>(u 2 + iv c 0)-€ 2 p 2 n } 2 
p m (co 2 + iv c LJ-{ 2 p 2 m ,)aL 



V c me ™fi cos ^ (JV C + 1/2) ] \ 
> SJk)S m/ (k) 1 p — — . 
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Near the fundamental-mode resonance a> m lp\ the ampli- 
tude a\ dominates and we can use single-mode approximation 
neglecting all other amplitudes. This leads to a simple result 



with 
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cos[ qi (N m - 1/2)] - (1 + JiOcos [qi (N m + 1/2)] 
giJn 



qi sin [qiN m ] - J u cos [q x (N m + 1/2)] 
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. With this result, one can ob- 



and q\ ~ y p- 

tain the oscillating phases and fields in the mesa. Also us- 



ing Eq. (10 1 and keeping only m = 1 in the sum, we ob- 
tain the coefficients which determine the oscillating mag- 
netic field inside the crystal, Eq. (|8j. For an isolated mesa 
on bulk crystal corresponding to the limit a,N c — > +oo, 
the following approximate result can be derived 38 JTn ~ 
VI — iVabO) (0.57 + 0.310 /€, suggesting the following pre- 
sentation Jii = VT — iv a bCL>Pi/£, where f}\ is the complex 
function of the order unity. The amplitude of the oscillating 
magnetic field on the top of the mesa can be represented as 



ipigi cos (pix) 



lo 2 - £ 2 p\ + iv c a> + Mi(a>) ' 
where the complex function J[\(<x>) 

(oj 2 - £ 2 p\ + iv c u))Pi 



(14) 



<*» ilM ^,(1 -cos [ qi (N m - 1)]) 

determines the damping of the resonance and its frequency 
shift in all quantities. 

Figure|3]illustrates the Josephson-frequency dependence of 
the oscillating magnetic field amplitude inside the mesa near 
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the cavity -resonance frequency for isolated mesa and for mesa 
array with a = 2L. One can see that for used parameters the 
corrections are weak. Above the main peak one can see a 
small dip caused to excitation of the almost uniform standing 
wave inside the crystal. To verify this, we also show the plot 
of the oscillating magnetic field at the bottom of the base crys- 
tal. The dip in the mesa field corresponds to the rather sharp 
peak of the field at the bottom of the crystal. Note also that 
the resonance is displaced to lower frequency with respect to 
the uniform cavity mode £n/L because plasma oscillations ex- 
cited inside the mesa are not uniform in the c direction. The 
width of resonance is mostly determined my leak of radiation 
into the base crystal. 

Figure [4] shows evolution of the resonance shape with vari- 
ation of the array period a. We can see that with increasing a 
the dip moves to smaller frequencies. The dip has maximum 
amplitude and located at the peak for a = 2.05L. It is inter- 
esting to note that the resonance in mesa is actually strongest 
when the dip is located below the peak at a = 2.1L. The rea- 
son is that in this case the wave excited at the peak frequency 
and the wave vector k - 2n/a is in decaying range, see Fig. 
[2] and, as consequence, the mesas loose less energy to radia- 
tion at the resonance frequency. Nevertheless, we expect the 
strongest interaction between the mesas and optimal condi- 
tions for synchronization when resonances coincide. 



IV. NUMERICAL SIMULATIONS 

To find the condition for synchronization between mesas 
and checked the stability of the synchronized state, we solve 
Eqs. ([T]i and (|2]i numerically for two mesas and numerical 
details are presented in Ref. |34| The number of junctions 
in the base is N c - 200 and in the mesa is N m - 50. We take 
v c = 0.02, v a b = 0.2 and £ = 266.5. To ensure that the resulting 
state is stable, we add an artificial weak white noise current in 
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FIG. 4. (color online) The frequency dependences of the oscillating 
magnetic field at the mesa top for different periods of the mesa array. 
Other parameters are the same as in the previous plot. 



FIG. 5. (color online) Phase coherence of plasma oscillations be- 
tween different stacks for different L m . The phase coherence is mea- 
sured by |gjx/«l defined in Eq. |T6|>, 
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FIG. 6. (color online) Snapshots of the electric field (first row), magnetic field (second row) and Josephson current sin(0„) (third row) in the 
whole system, (a) and (b) are obtained at the first and second cavity resonance for L,„ = L. (c) and (d) are results near the first cavity resonance 
for L m t nL. 



Eq. jlj in simulations, (l„(x, t)I„,(x', f')) = 10~ 5 S(x - x')6(t - 
t')6(n - n'). To study the coherence between different stack, 
we introduce an order parameter at edges of mesa 



N m 



1 



(15) 



where 0^ iL ^ R is the phase difference of n-th layer at the left 
(L) or right (R) edge of the fi-th mesa. The time average of 

lfy,L/«l> r n,LiR — Jo \ r ft,L/R\dt/T measures the phase coherence 
at the edges of the mesa. For coherent oscillations of phase 
difference fyx/fl = 1 an d f° r completely random oscillations 
r ti,L/R ~* when N m — > oo. To quantify the phase coherence 
between different mesas, we introduce a correlation function 



L/Rdt, 



(16) 



mode of the mesa a>\ — fn/L and the the peak at the right side 
corresponds to the second cavity mode a>2 = 2£n/L. When 
the frequency of the plasma oscillations in the mesa matches 
the cavity frequency a> m - mCn/L, \gj,L/R\ increases indicat- 
ing a tendency of synchronization between two stacks. When 
L m = nL, the phase coherence between different mesas be- 
comes maximal. This becomes clearer for the second cav- 



0.6 



where we have taken the left edge of the first mesa as refer- ^ 
ence. Similarly \gn,LiR\ measures the coherence between the 
phase at left or right edges of the fi-th mesa and the phase at 
the left edge of the first mesa, and the phase of gfi,L/R repre- 
sents the phase shift between them. 

Let us first consider two identical mesas with width L = 
0.3A C and with a separation L,„. They locate at the position 
L away from the edges of the base crystal. They are biased 
by the same current, I\ - h - hxu as shown in Fig. |6|a). 
The results for different L m is shown in Fig. [5] . The main 
peak at the left side corresponds to the fundamental cavity 
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FIG. 7. (color online) Voltage difference between two mesas AV = 
V2 - V\ when the mesas are biased by different current h = 1\+ A/ cxt . 
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ity mode, where two stacks do not synchronize at all when 
L m + nL. For L„, = 0.7 A c , the peak at 7 ext = 2.0I C is due to 
the cavity resonance inside the mesa. However the resonance 
occurs at smaller 7 e xt compared to that with L m = nL. The 
downshift is due to the strong radiation damping through the 
base crystal as shown in Fig. [3] 

The reasons for the better coherence when L m - nL are as 
follows. For the plasma oscillations uniform along the c axis 
qt = 0, the in-plane dissipation is absent and the plasma is 
damped by the weak dissipation along the c-axis according to 
Eq. However, for nonuniform oscillations with a finite 
wavenumber q^, the in-plane dissipation becomes dominant 
for N m ~ 10 3 , and the nonuniform plasma oscillations in the 
base crystal decays quickly. Therefore the interaction between 
two mesas is weak and the synchronization becomes difficult. 
When L m = nL, uniform cavity modes — are possible 
as shown in Fig. [6] (a) and (b). Two mesas then are strongly 
coupled through the base crystal and they are synchronized. 
When L m + nL, only nonuniform modes can be excited and 
the plasma oscillations in the base crystal is strongly damped 
by the in-plane dissipation. The amplitude of plasma oscilla- 
tions is small compared to that when L m - nL, see Fig. |6]c) 
and (d). Thus synchronization between mesas is hard to attain. 
Therefore the maximal synchronization is achieved when the 
position and size of mesas are commensurate with the stand- 
ing wave in the base crystal, because the nonuniform plasma 
oscillations decay quickly in the base crystal. 

Let us consider the phase shift of the gauge invariant phase 
difference between edges of mesas. As shown in Fig. |6ja), the 
supercurrent changes sign from the left edge to the right edge 
in the same mesa. This indicates there is a n phase jump at the 
center of the m esa, and the n phase kink is excited at the cavity 
resonanc d 32 l 33 l The n phase kink helps to pump energy into 
the plasma oscillations and the amplitude of the oscillations is 
enhanced sharply at the cavity resonances, as described by the 
term at the right-hand side of Eq. p). In Fig. |6|a) and (b), the 
plasma oscillations at the left/right edges have the same phase 
between different stacks when L m = (2n + 1)L. For L m = 2nL, 
there is n phase shift between two mesas to match the standing 



wave in the base crystal. 

In Fig. [7] the voltage of mesas with L m = L when they are 
biased by different current h = h + A/i ext is shown. At the cav- 
ity resonance when two mesas are synchronized, they have the 
same voltage despite that they are biased by slightly different 
current. Away from the cavity resonance, two mesas decouple 
from each other and they oscillate at different frequency. 



V. CONCLUSIONS 

We have investigated the synchronization of mesa array 
through the plasma oscillations in the base crystal. If one re- 
gards the mesa arrays and base crystal as a whole, the plasma 
oscillations inside depends on the configuration of mesa ar- 
rays as a result of geometrical resonance. The amplitude of 
the plasma oscillations and the tendency of synchronization 
is governed by the dissipation of the whole system, hence is 
determined by the configuration of mesa array. When the pe- 
riod of mesa array a close to the multiple integer of the mesa 
width L, a as nL, the dissipation is minimized and mesas are 
synchronized at cavity resonances of the whole system. Al- 
ternatively, one may treat mesa and base crystal separately. 
When the cavity resonance of mesa matches of that in the base 
crystal, mesa array is synchronized. Otherwise, the cavity res- 
onance of the mesas is suppressed by the strong dissipation 
due to the radiation into the base crystal, and the mesa array 
is not synchronized. Therefore the optimal configuration for 
synchronization is a « nL. The above picture is corroborated 
by both analytical calculations and numerical simulations. 
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